2. Model and simulation methodsWe consider a 2D system where
self-propelled particles (treated as disks) are enclosed by a soft boundary. Each particle has mass m, and diameter σ. In addition to the thermal Brownian motion at temperature T, each particle is propelled by a force with the amplitude F through its center. The orientation of the force,
rotates with time due to fluctuations. The soft boundary is mimicked by a ring of spring-connected passive beads (
) which are set to be identical to the active particles except being passive, i.e., without propelling force on them. All non-bonded interactions between particles and beads are given by the isotropic and purely repulsive Weeks–Chandler–Andersen (WCA) potential,[30]
with a cutoff at
, beyond which
, ε is the interaction strength. The bonded harmonic potential between boundary beads is
. We set the spring constant
and equilibrium bond length
to prevent the active particles from crossing the boundary. In the natural circular shape, the perimeter of the boundary
, the radius
and the enclosed area
.
The motions of particles and beads are described by the Langevin equations as follows:
Equation (
1) governs the translational motion. For the boundary beads,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn13.gif)
, the potential
Ui is composed of the non-bonded WCA potentials and bonded harmonic potentials, and only equation (
1) is required. In contrast, for the active particles,
Ui contains only the non-bonded WCA potentials and equation (
2) which depicts the coupled rotational kinetics of the propelling direction is necessary. The translational friction coefficient
![](cpb_26_10_100502/cpb_26_10_100502_ieqn14.gif)
with
D0 being the translational diffusion constant. The rotational diffusion rate
![](cpb_26_10_100502/cpb_26_10_100502_ieqn15.gif)
when only thermal noise is present.
![](cpb_26_10_100502/cpb_26_10_100502_ieqn16.gif)
and
![](cpb_26_10_100502/cpb_26_10_100502_ieqn17.gif)
are unit-variance Gaussian white noises.
We use the LAMMPS software to perform the simulations.[31] A square box of
with a periodic condition in both x and y directions is adopted. Reduced units are used by setting
,
, and
.
is the corresponding unit time. We set ε = 10 and γ = 10, which is large enough so that the motions of particles and beads are effectively overdamped. We treat
as an independent parameter because in many cases the rotational noise is athermal.[12,15] To highlight the influence of the particle activity, we set
to be a small value of
, inverse of which quantifies the persistence time of the active motion. For each case, it runs in a minimum time of
in time steps of
.
3. ResultsWe consider three initial area fractions of confined SPPs, ϕ = 0.05 (low), 0.2 (medium), and 0.4 (high), where
. Figure 1 shows the variations of the typical instantaneous morphologies with the increase of the propelling force on the particles for ϕ = 0.05 and 0.4. Apparently, in both cases, there are two regimes in the variation of the morphology with increasing propelling force. i) For small propelling force, the distribution of particles inside the closed boundary changes from homogeneous distribution to radially inhomogeneous aggregation, i.e., a concentration hole gradually forms in the center due to enhanced aggregation of particles near the boundary; the shape of the boundary roughly keeps circular; and apparently the local fluctuation of the boundary is increasingly suppressed. ii) For large active force, the distribution of particles becomes both radially and angularly inhomogeneous; the boundary is gradually elongated and deviates from a circle; this global shape deformation is the consequence of the inhomogeneous distribution of the particles along the boundary and the trend of their accumulating around the places of high curvature,[18] which forms positive feedback in the deformation of the boundary. In the following, we quantitatively analyze the variations of the size and shape of the boundary, the inhomogeneous distribution of particles, and their correlations.
First, we calculate the mean radius of gyration,
, of the boundary, based on the positions of its constituent beads, as a measure of its size for a ring polymer (Fig. 2(a)). In two dimensions, the radius of gyration is defined as
, where
and
are the two eigenvalues of the gyration tensor,
where
![](cpb_26_10_100502/cpb_26_10_100502_ieqn33.gif)
is the
α Cartesian coordinate of the
i-th bead on the boundary and
![](cpb_26_10_100502/cpb_26_10_100502_ieqn34.gif)
is the corresponding center-of-mass coordinate. At large driving forces (
![](cpb_26_10_100502/cpb_26_10_100502_ieqn35.gif)
, good linear relationships are found between
![](cpb_26_10_100502/cpb_26_10_100502_ieqn36.gif)
and
F for all the three initial area fractions of particles that we have calculated. The slope is larger for higher
ϕ. Figure
2(b) shows the variations of
![](cpb_26_10_100502/cpb_26_10_100502_ieqn37.gif)
with
ϕ under several (large) forces. The curves only slightly deviate from linearity. These nearly linear behaviors can be explained by the following simple scaling argument. Ignoring the noise for translational motion in Eq. (
1), the persistence length of a free active particle
![](cpb_26_10_100502/cpb_26_10_100502_ieqn38.gif)
. The ratio
![](cpb_26_10_100502/cpb_26_10_100502_ieqn39.gif)
can, to some extent, account for the degree of confinement.
[18] Figure
2(c) shows a dramatic decrease of the ratio at a small force, implying a steep transition from weak to strong confinement with the increase of force. The strong confinement causes the formation of the concentration hole in the center, i.e., the bulk particle density approaches zero (Fig.
1).
[18] At a large driving force, the thermal forces are negligible and it corresponds to the situation of strong confinement,
[12,18] i.e., all the active particles aggregate near the boundary and exert a force
F on it. Therefore, the total force,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn40.gif)
, exerted by the particles on the boundary, which causes its expansion, is proportional to the strength of the driving force and the number of confined particles. Ignoring the influences of inhomogeneous distribution of particles and the non-circular shape of the boundary, we have roughly,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn41.gif)
. This force is balanced by the tension of the boundary, i.e.,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn42.gif)
. Therefore, we have
![](cpb_26_10_100502/cpb_26_10_100502_ieqn43.gif)
, i.e.,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn44.gif)
is a linear function of
F or
ϕ.
At small F (inset of Fig. 2(a)), we find a crossover from rapid to slow expansion of the boundary with the propelling force around
. The initial rapid expansion is attributed to the aggregation of particles to the boundary, which effectively increases the average number of particles that exert a force on the boundary. We calculate the probability distributions of
for ϕ = 0.05 under different forces (Fig. 3(a)). The curve of F = 1 shows the narrowest distribution indicating the strong suppression of the local fluctuation of the boundary. This is further proved by the calculation of the relative fluctuation,
in Fig.
3(b) where minimum fluctuations for all particle area fractions are around
F = 1. The probability distribution at
![](cpb_26_10_100502/cpb_26_10_100502_ieqn48.gif)
is left-right symmetric. In contrast, the distribution at large force is asymmetric. These are the manifestations of the characteristics of the two regimes: radial expansion and suppression of local fluctuation of the boundary in regime one and the global shape deformation in regime two.
To quantify the global shape deformation of the boundary, we calculate the asphericity (deviation from a circle)[32,33] in Fig. 4. The asphericity in two dimensions is defined as
where
A is 0 if the shape is circular or isotropic. Figure
4(a) shows that
![](cpb_26_10_100502/cpb_26_10_100502_ieqn50.gif)
increases dramatically when
![](cpb_26_10_100502/cpb_26_10_100502_ieqn51.gif)
in the case of
ϕ = 0.05, while a slow-to-rapid increase is found for
ϕ = 0.2 and 0.4. This manifests the different manner (sharp transition at low
ϕ versus gradual transition at high
ϕ) in the crossover from regime one to regime two with the increase of propelling force. At low
ϕ, SPPs form a thin and sparse layer around the boundary at the end of regime one. Upon further increase of the propelling force, the SPPs easily slide along the boundary, which leads to a significant angular anisotropic distribution of particles and hence the shape deformation of the boundary. Therefore, the transition from regime one to regime two at low particle concentration is sharp. Presumably, the crossover happens around the condition satisfying the strong confinement, i.e.,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn52.gif)
which is consistent with the observed transition point
![](cpb_26_10_100502/cpb_26_10_100502_ieqn53.gif)
. In contrast, at high
ϕ, SPPs form a thick layer around the boundary and enhance the homogeneity of the angular distribution of particles. Consequently, in a large range of propelling force, the particles could not individually slide significantly but collectively along the boundary. Hence, the shape deformation keeps small even though the SPPs are strongly confined. Note that
![](cpb_26_10_100502/cpb_26_10_100502_ieqn54.gif)
for
ϕ = 0.05 shows a clear transition from the rapid to slow increase at a force around
F = 20. Examining the trajectories, we find the aggregation layer of particles along the boundary becomes thin and sparse for
ϕ = 0.05 and
![](cpb_26_10_100502/cpb_26_10_100502_ieqn55.gif)
(Fig.
1), and there is no appreciable adjustment of the thickness of the particle aggregation layers along the boundary upon the further increase of the force, which, otherwise, would contribute to the increase of
![](cpb_26_10_100502/cpb_26_10_100502_ieqn56.gif)
. Therefore, the slow increase at
![](cpb_26_10_100502/cpb_26_10_100502_ieqn57.gif)
is induced by the further expansion of the boundary without the contribution or enhancement from the further inhomogeneous angular distribution of active particles.
Time correlation function of asphericity
is a way to quantify the varying rate of the global shape (Figs. 4(b) and 4(c)). We find that the force-dependence of the decay is different for the low and high initial particle concentrations. For ϕ = 0.05, the
decays faster (i.e., the shape varies faster) at a larger force (Fig. 4(b)). The change of the decaying rate is appreciable. Taking the time at which
as the characteristic correlation time, we have that the correlation time changes from
for F = 5 to
for F = 80. On the contrary, the decay of
becomes slower upon the increase of force for ϕ = 0.4 (Fig. 4(c)), but apparently the change is small. The correlation time changes from
for F = 5 to
for F = 80. There are two competing factors in determining the rate of the shape variation: expansion of the boundary and the activity of the particles, both of which are enhanced with the increase of active force. The large size of the boundary slows down the shape variation while the large activity of particles leads to the rapid relocation and hence the fast shape variation. The expansion of the boundary with active force is weak for ϕ = 0.05 (Fig. 2(a)) and the enhanced activity of particles dominates and leads to faster decay of
at a larger driving force. In contrast, the expansion of the boundary is significant for ϕ = 0.4 and causes a much slower variation of the shape at a larger force. Notice that the correlation time for ϕ = 0.05 is larger than that for ϕ = 0.4 when
. We examine the trajectories and find that the possible reason is the different mechanisms for the particle relocation. For ϕ = 0.05 the layer of particles at the boundary is thin and the particles are relocated mostly through the motion along the boundary, while the layer of particles at the boundary is thick for ϕ = 0.4 and besides the motion along the boundary many particles near the center move across the central hole, which can definitely speed up the change of the angular distribution of particles.
To investigate the radial and angular inhomogeneous distributions of the enclosed particles and their coupling with the deformation of the boundary, we calculate the Gini coefficients[12]
and
by partitioning the space in radial and angular directions, respectively (see the inset of Fig.
5(a)), where
![](cpb_26_10_100502/cpb_26_10_100502_ieqn69.gif)
is the mean density,
ρi the number density of particles in the
i-th partitioned space, and
N the total number of partitioned spaces. For the radial Gini coefficient,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn70.gif)
, we choose the average position of the boundary beads as the center and divide the center-to-bead lines equally to form inner concentric boundaries which partition the space radially (
N = 40). We use the angular Gini coefficient,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn71.gif)
, to capture the inhomogeneous distribution of the particles along the boundary. To this end, we choose
N = 40 beads on the boundary with an interval of around 30 beads and draw circles with the chosen beads as the centers; the radius is empirically set to be
![](cpb_26_10_100502/cpb_26_10_100502_ieqn72.gif)
, where
L is the perimeter of the boundary. We also try other choices of
N and
![](cpb_26_10_100502/cpb_26_10_100502_ieqn73.gif)
, which only quantitatively changes the results. Figure
5(a) shows that for all initial area fractions,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn74.gif)
increases significantly when the force
![](cpb_26_10_100502/cpb_26_10_100502_ieqn75.gif)
and then gradually approaches a plateau. The turning point from regime one to regime two is around
![](cpb_26_10_100502/cpb_26_10_100502_ieqn76.gif)
or 2, above which most of the particles aggregate near the boundary. This is in regime one. Below the turning point, we find the angular Gini coefficient
![](cpb_26_10_100502/cpb_26_10_100502_ieqn77.gif)
decreases (Fig.
5(b)). We ascribe this decrease to the suppression of fluctuations of both the particles and the boundary when particles aggregate near the boundary. At a large
F, all values of
![](cpb_26_10_100502/cpb_26_10_100502_ieqn78.gif)
increase with force, manifesting the enhanced angular inhomogeneous distribution of particles, which corresponds to the regime two. The crossovers from regime one to regime two are found to be different at low and high values of
ϕ. For
![](cpb_26_10_100502/cpb_26_10_100502_ieqn79.gif)
increases immediately when
![](cpb_26_10_100502/cpb_26_10_100502_ieqn80.gif)
, i.e., the crossover is sharp. In contrast,
![](cpb_26_10_100502/cpb_26_10_100502_ieqn81.gif)
keeps small and does not rise appreciably until
![](cpb_26_10_100502/cpb_26_10_100502_ieqn82.gif)
for
ϕ = 0.2 and
![](cpb_26_10_100502/cpb_26_10_100502_ieqn83.gif)
for
ϕ = 0.4, i.e., the crossover is gradual. These behaviors are in consistent with the variations of asphericity in Fig.
4, manifesting a strong correlation between the shape deformation of the boundary and the inhomogeneous angular distribution of particles.
Finally, we calculate the local curvature κ of the boundary and also the local pressure on it by statistics. The local curvature
, where (
) is the coordinate of a boundary point, and
and
are the second and first derivatives with respect to x. The local curvatures are calculated at the bead positions along the boundary in the discrete form of the derivatives. The local pressure exerted by active particles is evaluated by the local force along the local normal direction
of the boundary divided by the local arc length. The statistical data for ϕ = 0.05 are shown in Fig. 6(a). It is found that the local pressure is almost independent of the local curvature at small forces (
), which corresponds to the small asphericity of the boundary as a whole (Fig. 4(a)). For
, local pressure increases with the increase of the local curvature, which accords with the observation on fixed boundary.[15,18] The correlation between local pressure and local curvature is directly visualized in Fig. 6(b), which shows the heat map of the local pressure and curvature along the contour of the boundary for the case of F = 80 and ϕ = 0.05.